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ABSTRACT 

In this publication the construction of an automatic algorithm to subtract infrared divergences in 
real QCD corrections through the Catani-Seymour dipole subtraction method [lj is reported. The 
resulting computer code has been implemented in the matrix element generator AMEGIC++ [2j. This 
will allow for the automatic generation of dipole subtraction terms and their integrals over the one- 
parton emission phase space for any given process. If the virtual matrix element is provided as well, 
this then directly leads to an NLO QCD parton level event generator. 



tanjuOtheory . phy . tu-dresden . de 
2 f rank . krauss@durham .ac.uk 



Contents 



1 Introduction 

2 Brief review of the Catani-Seymour formalism 

2.1 NLO cross sections and the subtraction procedure 

2.2 Generalisation to hadronic initial states 

2.3 Observable-independent formulation of the subtraction method 

2.4 The dipole subtraction functions 

2.4.1 Final state emitters with final state spectators 

2.4.2 Final state emitters with initial state spectators 

2.4.3 Initial state emitters with final state spectators 

2.4.4 Initial state emitters with initial state spectators . . . . 

2.5 Integrated dipole terms 

2.5.1 Phase space factorisation 

2.5.2 Full result 

2.6 Freedom in the definition of dipole terms 

3 Implementation in AMEGIC++ 

3.1 ME generation 

3.1.1 Amplitudes 

3.1.2 Process structure and organisation 

3.2 Generation of CS dipole terms 

3.2.1 Colour and spin correlations 

3.2.2 Organisation and process management 

3.3 Generation of the finite part of integrated dipole terms 

3.3.1 Analytical structure of the full result 

3.3.2 Implementation and Organisation 

3.4 Phase space integration 

3.4.1 Importance sampling and Multi-channel integration . . 

3.5 Cuts and analysis framework for NLO calculations 

4 Checks of the implementation 

4.1 Explicit comparisons 

4.2 Test of convergence for the real ME 

4.3 Consistency checks with free parameters 

5 First physical applications 

5.1 Three-jet observables at LEP 

5.2 DIS: e'p -> e~ + jet 

5.3 W~ production at Tevatron 



6 Conclusions and outlook 

7 Acknowledgments 

A Insertion operators etc. 



1 Introduction 



Perturbative calculations form one of the best understood methods to provide predictions for the 
behavior of a Quantum Field Theory and to compare them with experimental results. Many of the 
methods applied in such calculations have found their way into textbooks already decades ago, see 
e.g. [3]- [Jj. Typically, the perturbation parameter is related to the coupling constant of the theory 
in question, which in most cases indeed is a small quantity. This also implies that the corresponding 
fields may asymptotically appear as free fields and thus are the relevant objects of perturbation 
theory. Obviously, this is not true for the strong interactions, i.e. QCD, where the fields, quarks 
and gluons, asymptotically are confined in bound states only. This is due to the scaling behavior of 
the coupling constant of QCD, as, which becomes small only for large momentum transfers, see for 
instance |8)9] . It is the confinement property that to some extent restricts the validity of perturbative 
calculations in QCD to the realm of processes characterized by large momentum transfers or by other 
large scales dominating the process, such that factorization theorems can be applied [9]- [TT] . 

Typically, for most of the relevant observables in particle phenomenology, the leading term of the 
perturbative expansion can be related to tree-level diagrams. In the past years, the calculation of 
these terms has been fully automated and a number of tools capable of dealing with up to eight to 
ten external particles without any significant user interference have emerged |2|12j - |16| . However, 
for many practical purposes, tree-level calculations are not sufficient. This is due to a number of 
reasons: first of all, many measurements aim at the extraction of fundamental parameters. However, 
in Quantum Field Theories, parameters are subject to corrections, which usually exhibit ultraviolet 
divergences. These divergences are dealt with through the renormalization procedure, which can 
be done in a scheme- and scale-dependent way only, see e.g. [6 1 1 7| . Therefore, in order to extract 
parameters from the comparison of a (perturbative) calculation with experimental data, the calcu- 
lation itself must contain the same kind of quantum corrections necessitating their renormalization. 
Second, it should be stressed that in tree-level calculations, there are some choices to be made, con- 
cerning the scale at which inputs such as the coupling constant, quark masses or parton distribution 
functions are taken. In principle, different scale choices are equivalent, and renormalization group 
theory guarantees that, when taking into account all orders, the effect of scale choices vanishes. At 
leading order (LO), however, their impact may still be significant, such that tree-level calculations 
merely give the order of magnitude for corresponding cross sections etc.; a prime example for this 
is the production of a Higgs boson in gluon fusion processes, where only the next-to-next-to lead- 
ing order correction significantly reduces the scale dependence and produces a stable result |18|19| . 
Thus, aiming at any more precise prediction, higher-order calculations are a crucial ingredient of 
phenomenological analyzes. 

But although indispensable, so far there is no fully automated tool available for QCD calculations 
at next-to leading order (NLO), i.e. at the one-loop level. This is because a true NLO calculation 
is certainly much more complex than a leading order (LO) one. First of all, some of the essential 
ingredients, namely the loop or virtual contributions are not under full control yet. In general, up 
to now calculations of these corrections to physical processes are limited to contributions containing 
five- and in some cases six-point functions, see for example [20]- |25| . But even to reach the level 
of known scalar master integrals is far from being trivial; the tensor reduction necessary for this 
step |26| results in a proliferation of terms with non-trivial cancellations among them, which render 
the implementation in a computer code a major effort. On the other hand, some of the loop 
corrections exhibit not only ultraviolet divergences to be renormalized, but also infrared divergences. 
They also need to be regularized, but then they must be canceled against similar infrared divergences 
stemming from the real contributions. This basically translates into canceling divergences in phase 
space volumes of different dimensionality. The cancellation in fact is one of the most important 
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consequences of the Kinoshita-Lee-Nauenberg or mass factorization theorems [27J- [29J. However, in 
order to practically achieve the cancellation, the real infrared divergences also need to be regularized. 
Essentially, there are two ways of doing this. 

One method, also known as phase-space slicing [30]- [35], bases on dividing the phase space of the 
additional real emission into an infrared-safe (hard) and a infrared-divergent (soft) region. The 
division is usually performed by subjecting pairs of particles to an invariant mass criterion. Then, 
the soft region is integrated analytically in d dimensions. Typically, in this step, the helicity- 
summed matrix element squared is approximated by its double-pole (or eikonal) limit. The result 
of the analytical integration will contain single or double poles of the form l/(d — 4) or \/{d — 4) 2 , 
respectively. They typically are accompanied with logarithms of the invariant mass criterion. Such 
logarithms, but with opposite sign, also appear in the numerical evaluation of the full matrix element 
squared for real emission in the hard region of phase space, performed in 4 dimensions. In principle, 
these two potentially large contributions (logarithms of a potentially small quantity) originate from 
the unphysical division of the phase space and should thus cancel. Therefore, the key issue thus in 
phase space slicing is to adjust the parameters of the procedure such that the dependence on the 
slicing parameter is minimized. So far, this adjustment has been done manually only and this is one 
of the reasons why other methods have become more popular with practitioners of NLO calculations. 

Such alternative methods of dealing with the real infrared divergences base on directly subtracting 
them [l|36j - [32J3. At NLO level the subtraction in all methods is performed such that the additional 
particle is added to the leading order matrix element in a well-defined way through terms which, on 
one hand, exhibit the correct divergent behavior in the soft and collinear limit, and, on the other 
hand, can easily be integrated over the full d-dimensional phase space of the extra particle. The 
idea is then that the so subtracted matrix element squared is finite and thus can safely be integrated 
numerically in 4 dimensions. On the other hand, the subtraction term is added to the virtual bit 
and integrated analytically in d dimensions. Again, it exhibits single or double poles of the form 
l/(d — 4) or l/(d — 4) 2 , respectively. These poles again cancel the infrared poles of the virtual 
contributions. The fact that there are universal subtraction terms, i.e. terms which will cancel the 
infrared divergences in a process- independent manner, is one of the main reasons why subtraction 
methods have become increasingly popular in past years and why they have been used for many of 
the state-of-the-art calculation of NLO corrections to physical processes, like for instance [20]- [24| . 

The universality of the subtraction terms also allow for an automated treatment of real infrared di- 
vergences. It is the subject of this publication to report on a fully automated, process-independent 
implementation of one of the popular subtraction procedures, ready for use in realistic NLO calcu- 
lations. Therefore, the outline of this paper is as follows: In section 2, the anatomy of QCD NLO 
calculations will be formalized in a more mathematical language and the chosen subtraction method, 
the Catani-Seymour dipole subtraction [1] will briefly be reviewed in its original form for massless 
particles. Although its extension to massive particles [10J is straightforward from an algorithmic 
point of view, this paper concentrates on the massless case only. In section 3, the fully automated 
implementation of the corresponding massless dipole subtraction of arbitrary matrix elements into 
the matrix element generator AMEGIC++ [2] will be presented in some detail. Some simple tests 
of the implementation will be discussed in section 4, before some physical applications and the 
comparison with results from the literature will round of the presentation in section 5. 

1 Subtraction methods for the NNLO case have been presented in [43]- [51] . 
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2 Brief review of the C at ani- Seymour formalism 



2.1 NLO cross sections and the subtraction procedure 

Cross sections at NLO precision are given by 

a — a + a , [L) 

where the LO part a LO is obtained by integrating the exclusive cross section in Born approximation 
over the available phase space of the m final state particles and, eventually, over the Bjorken-x of 
incident partons. Ignoring this additional complication for the sake of a compact notation, The LO 
cross section is thus given by 



a LO ^ / j(4)_B 



d< 4 V B , (2) 



where 



d^a B = d^^\M m \ 2 Fp. (3) 

Here, d( 4 )<E>( m ) denotes the phase space element of m particles, taken in four dimensions, A4 m is 
the matrix element for the process under consideration, and Fj 11 " 1 is a function of cuts defining the 
jets etc.. As already indicated, here and in the following, the superscripts in the integral denote 
the dimensionality of the integration. In order to obtain a meaningful result to be compared with 
experimental data, typically isolation cuts are applied on the outgoing particles, which may also serve 
the purpose of keeping the integral finite. A typical criterion for example is to identify outgoing 
partons with jets and thus apply jet definition cuts on the partons such that they are all well 
separated in phase space. Anyway, the cuts will not be stated explicitly in the integral, but they are 
understood implicitly with the integration, including suitable generalisations in d dimensions, where 
necessary. Thus, the integration of the Born level cross section can directly be carried out in four 
space-time dimensions, as indicated in the equation. 

In view of the dipole subtraction formulae, it is useful to introduce at this point bras and kets 
m (l, . . . , m'\ and |1, . . . , m') m . They denote states of m final state partons partons labeled by 1 to 
m! and are vectors in colour and helicity space. Introducing, in a similar fashion, vectors for the 
spins and colours, matrix elements thus can be written as 

Mm" 1 = (m(ci,...,C m |<g) m (si,...,S m |)|l,...,m) m . (4) 

Therefore, in this notation, the matrix element squared, summed over final state colours and spins 
reads 

\M m \ 2 = m (l,...,m\l,...,m) m . (5) 

The NLO part of the cross section consists of two contributions, each of which increases the order 
of as- First, there are emissions of an additional parton, i.e. real corrections, denoted by d^ rf V R . 
Second, there are virtual (one-loop) corrections to the born matrix element, here denoted by d( d V v . 
Thus, 

^nlo = r d ( d)a NLo = r d ( d ) a R + r d( v (6) 

J J m+l Jm 
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The two integrals on the right-hand side of Eq. (|6|) are separately infrared divergent in four dimen- 
sions, and are therefore taken in d dimensions. For the real correction, the divergences arise when the 
additional parton becomes soft or collinear w.r.t. some other parton, leading to on-shell propagators 
in the matrix element. For the virtual correction, the divergence comes with the integration over 
the unrestricted loop momentum, such that again a propagator goes on-shell. As already stated 
in the introduction, now the celebrated theorem of Kinoshita, Lee and Nauenberg |28|29] comes to 
help and guarantees an exact cancellation of two divergent contributions, thus keeping their sum 
finite! Setting d = 4 + 2e in the following, the divergences will manifest themselves in double and 
single poles, i.e. as 1/e 2 and 1/e, respectively. In principle, cancellation of the poles then solves 
the problem; in practice, however, the direct applicability of the equations above to real physical 
processes is limited since analytical integration over a multi-particle phase space in d dimensions 
with cuts in many cases is beyond current abilities. 

Therefore, a detour has to be taken. The idea is to construct a subtraction term for the real emission 
contribution, which encodes all of its infrared divergences, but can analytically be integrated over 
in d dimensions. In this way the infrared pole structure of the real part with its 1/e and 1/e 2 
poles is exhibited and cancels the corresponding virtual contributions. Subtracting this term from 
the real emission contribution and adding it to the virtual corrections then eliminates the infrared 
divergences in both parts. The subtracted real matrix element squared then is finite and thus its 
full (m + l)-particle phase space can safely be integrated over in four dimensions. In this way, the 
subtraction term aims at an infrared regularisation of the two contributions at integrand level. 



a 



NLO 



R 



d< d V A + 



m+1 J m+1 

d ( 4 V R -d( 4 ><7 A 

m+1 L 



d (d)_A 



m+1 



+ 



+ 



m+1 



J m 



d(V 



V 



(7) 



The catch of the subtraction method now is that the subtraction terms can be obtained from the 
Born terms in a straightforward way and that only the phase space integral of the extra particle 
has to be taken in d dimensions, while the phase space for the remaining m particles can be taken 
in four dimensions. This is similar to the way, the loop terms are evaluated. There, only the loop 
integration is performed in d dimensions, whereas the phase space of the outgoing particles is done 
in four dimensions. Therefore, the final structure reads 



a 



NLO 



m+1 



d (4) 



a 



A 



v 



+ 



loop 



d («0 



a 



(8) 



Both integrands now are finite, allowing all integrations to be performed numerically. In contrast to 
some other regularisation methods (like, e.g., phase space slicing) the subtraction method does not 
rely on any approximation and does not introduces any ambiguous and/or unphysical cut-off scales 
etc., as long as the integration of d^V A can exactly and analytically be performed. 

In pQ a general expression for d^a A has been presented, called the dipole factorisation formula, 
allowing to write 



d (d) 



a 



X; d( v ® d( rf v dipolc 

dipoles 



(9) 



2 In fact, this is only guaranteed for infrared-safe quantities. More specifically, if defines jets in terms of the 
momenta of an m-parton final state (taken at Born level), infrared safety demands that F < ^ n+V> — > p( m ) m cases w here 
the m+1- and m-parton configurations become kinematically degenerate. 
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such that, symbolically, 



I m+1 



dipoles 



ipole 



(10) 



where 



dipole 



:n) 



dipolos 



Here the sum of the dipole terms Vdi po ie contains all soft and collinear divergences of the real emission 
pattern. This factorisation formula is suited for any process with massless partons, and fulfills all 
the requirements mentioned above. An extension to massive partons has been presented in [40] . 

However, as already mentioned in the introduction, in this publication only the massless case will be 
considered. In order to provide a self-contained description, all necessary analytic expressions will 
be listed in this publication. 



2.2 Generalisation to hadronic initial states 



The cross sections discussed so far were given for point-like initial states. For cross sections in hadron 
collisions, however, the differential cross sections above must be convoluted with parton distribution 
functions (PDFs): 



(12) 



a.b 



Here the subscripts on the cross section denote the flavours of the incoming partons; for the total 
cross section a sum over them has to be performed. For the NLO part, now the higher-order 
corrections residing in the PDFs must be taken care of. This is done by supplementing the NLO 
part with a collinear subtraction term daQ, such that 



^ b h °{pa,Pb,^F)= I d(<M(Pa,P&)+ / d^a v ab {p a , Pb ) + [ d^ag( Pa , Pb , /4) 

J m+1 J m Jm 

This new term contains collinear singularities, incorporated in 1/e-terms and reads 



(13) 



(zp a ,Zp b ) 



S bd S(l - z) ( -\ (^pj Pac { z ) +K F J-(z] 



+5 ac 5{l - z) 



1 



\14 



2\ e 



P bd (z) + Kll(z 



(14) 



The collinear subtraction term is factorisation-scale and scheme dependent. This scheme dependence 
resides in the terms K F , which, for the common MS-scheme vanish, i.e. in this scheme all terms 
j(F S. = o However, this scheme dependence cancels similar terms in the PDFs such that, taken 
together, the full hadronic cross section again is scheme-independent. 
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In the case of incoming hcicLrons, the subtraction method is ctpplied. to o - ^^ 1 ® (jj a ^ Pb,^ F ) as described 



before, with the only difference that in this case the singularities of dcK only cancel in the sum 



loop 



(15) 



<E = 



2.3 Observable-independent formulation of the subtraction method 

Up to now, the da denoted cross sections in a broad sense. To be a bit more specific consider the 
following expression for a cross section at Born-level and the corresponding next-to leading order 
expression: 



a 



NLO 



d^ m \ Pl ,..., Pm ) M( m \ Pl ,...,p m ) F^(pi 



,Pr> 



d^ m+1 \ Pl ,...,p m+1 ) M^ m+1 \ Pl 



■•>Pm+l, 



F( m+1 \ Pl ,..., Pm+1 ) 



+ J d^ m \ Pl ,...,p m ) \v^\ Pl ,..., Pm )\ 2 F^\ Pl ,...,p m ) 



(16) 



where d$( n ) represents an n-particle phase space element, and M^ m \ M( m+1 ) and are the 

LO matrix element, the NLO real matrix element and the NLO virtual correction matrix element, 
respectively. F^ is a function that defines a cross section or an observable in terms of the n- 
parton momentum configuration. In general, the function F may contain ^-functions (to define cuts 
and corresponding total cross sections), 5- functions (defining differential cross sections), kinematic 
factors or any combination of these. 

However arbitrary this sounds, there is a formal requirement on this function F, namely that in 
the soft and collinear limits, i.e. for cases where one parton becomes collinear w.r.t. another one or 
where one parton becomes soft, the function _F( m+1 ) reduces to F^: 

F( m+l \ Pl ,..., Pl = Xq,...,p m+1 ) -» F( m \ Pl ,...,p m+1 ) for A^O 

F (m+1) {Pi,-,Pi,-,Pj,-,Pm+i) F (m \pi,-,P,-,P m +i) for pi -> zp, pj -► (z - l)p 

F^(p u ...,p m ) (17) 

The first two conditions define infrared-safe observables - to phrase it intuitively this means that 
such infrared-safe quantities must not be altered by additional soft or collinear activity. The last 
condition above is required to properly define the Born cross section. 

Applying the subtraction method to the NLO-part of Eq. (Tl6l) results in 



a 



NLO 



(j$(m+l) 



M^ m+1 \ Pl , ...,p m+1 ) F^ +1 \ P1) ...,p m+1 ) 



+ J d$ (m) 



V ij,k{pi, ...,p m +l) F^ m \p 1 ,..,p ij ,p k ,..,p m+1 ) 

v^(pt,...,p m )\ 2 +(^jmv lhk {p 



1) — ,Pm+l, 



F( m \ Pl ,..., Pm ), 
(18) 
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where d[p] is the phase space element for the 1-parton phase space. 

In order to have an identity between the subtracted terms and the added term, both the (m + 1)- 
parton contribution and the m-parton contribution have to be subjected to the same function F. To 
be able to perform the integration over the one-parton phase space independent of the observable 
this function therefore must be F^ m \ In the case of the (m + l)-parton contribution _F( m ) is applied 
to the m-parton configuration, generated by corresponding mapping given in the prescription of the 
dipole function. 



2.4 The dipole subtraction functions 



The universality of the soft and collinear limits of QCD matrix elements are the basis for the 
construction of the dipole subtraction terms. In both limits any matrix element squared for m + T 
partons factorizes into an m-parton matrix element times a (singular) factor. 

To be specific, consider first the soft limit of the matrix element, given by the momentum pj of 
parton j becoming soft, i.e. = Xq^ with A — > []. Then, employing 



PiPk 



PiPk 



(piq) (Pkq) (piq) [(pi +Pk)q] 

the soft limit reads 

m+l(l, • • • ,J, • • • ,m + 1|1, 



+ 



PiPk 



[ipi +Pk)q]{pkq) ' 



(19) 



,J, 
1 



-87r^ e a s ^ 

i,k^i m 



,J, 
h. 



,m + 1 



PkPiTk ■ Tj 



{piq)[(pi +Pk)q] 



, m + 1 



.(20) 



In a similar way the limit where two partons i and j become collinear is defined through pj 
(1 — z)/z pi. In this limit the (m + 1) parton matrix element can be rewritten as 



1+1(1, ...,m + 1||1, ...,m + l) m +i 



1 



PiPj 



1, . . . ,m + 1 



P(ij),i{z,k ± ) 



1, . . . ,m + 1 



(21) 



where, again, the Pujy^z, k±) are the well-known Altarelli-Parisi splitting functions. 

Then, the actual dipole function generating the limit, where one of the partons i, j of a m + 1-parton 
configuration becomes soft or both partons become collinear to each other, symbolically has the 
following structure: 



ij,k 



n (l, ...,k, ...,m||l, ...,k, ...,m) m O V ijtk , 



(22) 



with the non-singular m-parton matrix element m (...\\...) m and the operator Vj^fe, describing the 
splitting of the parton (ij). Here, and in the following, the splitting kernels Vy fc are matrices in 
the helicity space of the emitter. The dipole function also involves a third parton as 'spectator'. 
This parton in fact is identical with the colour partner k in the soft limit, Eq. (I20p . The form of the 
subtraction means that kinematically, 3^2 mappings are considered 



Pi,Pj,Pk -> Pij,Pk 



(23) 



such that all involved partons are allowed to remain on their mass shells. 

In general the splitting parton (called 'emitter') and the spectator can be both, initial and final state 
particles. This discriminates four different types of dipole functions, displayed in Fig. [TJ 
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Figure 1: Classification of dipole functions 
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The full subtraction term for any matrix element with (m + 1) partons in the final state is given by 
the sum of all possible dipole functions. For the most general case with two partons in the initial 
state, therefore 



v a a + S + Yl vai ' b + ( a ~ b ) 

i^j k+i i 



(24) 



In the following the explicit expressions for the dipole functions will be listed. The corresponding 
one-parton phase space integrated subtraction terms are discussed in Sec. 



2.4.1 Final state emitters with final state spectators 



The dipole contribution "Dijk for the singular limit p%-pj — > 0, where all three involved partons are 
in the final state, is given by 



T^ij,kiPU ■ ■ ■ ,Pm+l) 
1 / 



2pi -pj 



,k,...,m+l 



Tfc • T,; 



rp2 

ij 



-V- • 

▼ 1.1 



ij,k 



, m + 1 



(25) 



It is obtained from an (m + l)-parton matrix element by replacing the partons i and j with a single 
parton ij, the emitter, and the parton k is replaced by k, the spectator. The flavours of emitter and 
spectator are assigned as follows: The spectator k remains unchanged, and the emitter ij is defined 
by the splitting process ij — ► i + j. The product of colour charges in the numerator of Eq. (I25p 
introduces an extra colour correlation in the m-parton matrix element. 

The kinematics of the splitting are described by the following variables 



Uij,k 



Pipj 



Pipk 



1 



(26) 



Pipj + Pjpk + Pkpi Pjpk + Pipk 

and to obtain the momenta ij and k in the m-parton configuration the following map is being used: 



-A' 
Pk 



1 



u 

-Pk 



P 



1 Vij,k ^ J 1 Vij,k 

Obviously, four-momentum conservation is exactly fulfilled, i.e. 

a ii ii ~u . ~u 
Pi+Pj +Pk=Pij +Pk 

and all partons remain on their mass shell, 

0. 



2 2 2 -2 
Pi = Pj =Pk= Pij 



Pk 



(27) 



(28) 



(29) 



The splitting matrices, which are related to the d-dimensional Altarelli-Parisi splitting functions, 
depend on the spin indices of the emitter parton. For the case of a quark splitting (using helicity 
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indices s and s') the kernel is a matrix in helicity space, whereas for gluon splittings (to a quark- 
anti-quark pair or to gluons), the splitting matrices are given by Lorentz tensors. This yields 



(s\Vq ig:j ,k(zi;yij,k)\s') = s-K^ 2e a s c F 

(v\V qi q jt k(zi]yij,k)W) = 8irn 2e a s T R 
(Ml^ w ,fc(ij;ytj,jfe)H = 16nfJ? € asCA 



1 - Zi(l - y ij)k ) 



(1 + zi) - e(l - Si) 



PiPj 



(ziPi - ZjPjY(ziPi - ZjPjY 



+ 



1 - Zi(l - y ijik ) 1 - zj{l - y ij7 k) 



+(1 - e) (^pi - ZjPjY(ziPi - ZjPj) 1 

PiPj 



(30) 



respectively. The dipole terms given in this section are sufficient for the subtraction procedure in 
the case of non-hadronic initial states such as e _ e + -annihilation. 



2.4.2 Final state emitters with initial state spectators 

For the case of an emitting final state parton, the presence of an initial state spectator results in 
additional contributions to the singular limit Pi-pj — > of the full m + 1-parton matrix element. The 
corresponding dipole terms in this case are given by 



Vfjipi, . . . ,p m+ i;Pa, ••) 



1 1 



2pi 'Pi Xi h a m ,a 



l,...,ij,...,m + l;a,.. 



T • T 

ij 



l,...,ij,...,m + l;a,.. 



The kinematic variables now read 
_ i PiPj s . 



PiPa 



1 — £j 



{Pi+Pj)Pa' PjPa+PiPa 

and the momenta of the m-parton configuration are obtained by the map 

Pa = Xij^aPa 1 , Pij = Pi + Pj ~ ~ 2ty, )j# ■ 

Again, four-momentum conservation is trivially fulfilled and the partons remain massless. 
The corresponding splitting functions used in Eq. (|3T|) read 



(31) 



(32) 



(33) 



{s\V"(zi;x ijta )\s') = Snn 2e a s C F 



1 2fj ~\~ (1 %ij,a) 



[1 + z^ - e(l - Zi) 



(li\V£q.(zi;x ijja )\v) = Sjrfjra s T R 



PiPj 



[ZiPi - ZjPjY(ZiPi - ZjPj) 1 



+ 



1 ~H (1 1 (1 ^ij,a 



+(1 - e) (%>j - ZjPjY{ziPi - ZjPj) 1 

PiPj 



(34) 



12 



2.4.3 Initial state emitters with final state spectators 



The next type of dipole function now covers initial state singularities p a -pi — > with final state 
spectators, given by 



■ ■ ■ ,p m +i;Pa, ■■ 
1 1 



^Pa'Pi %ik,a 



1, . . . , k, . . . ,m + 1; ai, .. 



Tz, • T„ 



- k - CLi -trgj 
rp2 V k 

ai 



1, . . . , k, . . . , m + 1; ai, .. 



m,a 

(35) 



The parton ai, which enters into the m-parton matrix element on the r.h.s. of Eq. (J35J) is given by 
the splitting of the initial state parton a — > ai + i. The relevant kinematic variables in this case are 



%ik,a 1 



PiPk 



U; 



PiPa 



l-U k 



{Pk +Pi)Pa ' " PiPa+PkPa 

and the momenta for the m-parton configuration are obtained by 

Pai = x ik,aP^ , Pk=P k +P' 1 i - (! - X ik ,a)P^ ■ 

The splitting matrices V^* in Eq. (|35|) are 

2 



(36) 



(37) 



{s\V^( Ul ;x ikA )\s') 
(s\V k 9 ^(ui;x tk , a )\s'} 
(»\V k qaqi ( Ui ;x ik>a )\v) 

(v\V k 9i9a (ui;x ikta )\v) 



1 x ik,a ~i~ ^i 



8717J e a s C F [1 - e - 2x ik>a (l - x ik>a )\ 5 ss i , 



WTr^asCA 



9^ •Eik,a ~i~ 



- 6(1 - 


Xik,a) 


fiss' , 




(pi_ _ 


»y 


' ' Pi_ _ 


Pk\ 




UkJ 




u k ) 



1 



1 + Xi ka (\ 2-ifc,a)^ 



+ (!-£) + 



1 a 

UjU k 1 - X ikj q fPj_ _ PkY (Pj_ _Pk_ 
PiPk Xik,a \Ui U k J \Ui U k 



(38) 



The three dipole types discussed up to now (FF, IF, FI) are sufficient to construct the subtraction 
term da A for processes with exactly one initial state parton, i.e. DIS configurations. 



2.4.4 Initial state emitters with initial state spectators 

The remaining dipole function, only required by processes with two initial state partons, covers the 
case where both, the emitter and the spectator, are initial state particles, 



V ai ' b ( Pl ,... ,p m+1 ;p a ,p b ) 



11 /- ~ ~ , 
(l,...,m + l;ai,b 

Z Pa-Pi x i,ab mnh \ 



1, . . . , m + 1; ai, b 



m,ab 



(39) 



To describe the splitting, the following kinematic variables are used 



_ _ PiPa + PiPb ~ _ PaPi 



PaPb 



PaPb 



(40) 
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The construction of the m-parton kinematics for this dipoles differs from the other three cases. The 
reason is that in this case the emitter and the spectator are fixed to remain along the beam axis. 
Therefore all final state momenta (not only momenta of QCD partons) are transformed according 
to the map 



{K + K) 2 



K 2 



where 



lf=p£ + p£-# and £"=p&+p£. 



(41) 



(42) 



The momentum of the spectator p b remains unchanged. The transformation above can also be 
interpreted as applying a rotation and a boost turning initial state momenta back to the beam axis 
after a mapping similar to the first three cases of dipole functions. Indeed it can be shown that the 
transformation of final state momenta in Eq. (I41|) is just a Lorentz transformation. 

However, in this case, the splitting matrices read 



(s\V^<\x ijah )\s') 
(s\V^ b (x hab )\s') 

{n\V^ b {v l]X ^ a )\v) 

^\V^ b {v i]Xikta )\u) 



1 Xi ab 



(1 + Xi tab ) - e(l - Xi. 



ah j 



8tt/j, e a s T R [1 - e - 2x i<ab (l - x i>ab )] 6 S 



87r/j? e asCF 



g^XLab + - 



2 1 %i,ab 



r div 



Vi Pi -Pb X iyab 
Xi ab 



(pi - vip k Y (pi - Vip k ) v 



1 Xi,ab 



~t~ Xi,ab(X X 



i,ab, 



+(1 - e)-^ — (Pi - ViPkY (Pi - ViPk) 1 

Vi pi -p b x iAb 



(43) 



2.5 Integrated dipole terms 
2.5.1 Phase space factorisation 

In order to combine the poles of the subtraction function and the virtual matrix element the sub- 
traction function has to be integrated analytically over the one-parton phase space of the respective 
splitting. The rules for the momentum mapping from 3 to 2 parton phase spaces have been con- 
structed in Sees. I2.4.H|2~01 such that the corresponding phase space exactly factorizes. 

As an example, and in order to fix the notation, the case of a final- final dipole, V^j k, will be discussed 
in the following. There, the three-particle phase space for the partons i, j and k (all other partons 
are not affected by the splitting and will be omitted) in d dimensions is given by 

d(/)(pi,Pj,p k ;Q) = ( 2 t)rf-i S+ ( p ^ (2txY- 1 (2ir)d-i 6 +(Pk)( 27T ) d6(d) (Q-Pi- Pj ~ Pk) ■ 

(44) 

This can be factorized in terms of the mapped momenta, such that 

d4>(pi,pj,p k ;Q) = d<f)(pij,p k ;Q)[dpi{pij,p k )] , (45) 
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where [dpi(pij,Pk)], written in terms of the kinematic variables defined in section [2.4.11 reads 



[dpi(pij,Pk] 



(2p ij p k ) 1 ^ drt d -V 

16vr 2 (2vr) 
izi(l-Zi))- e (l 



1-2<E 



dzi dy ijik 9(%i(l - Zi))0(y ijik (l - y ijik )) 



(46) 



Within the dipole function only the splitting function itself depends on the variables Zi and yij,fc- 
Thus, the integration in d dimensions can be performed once and for all, independent of the specific 
scattering process under consideration. The result of the integration for each splitting type can be 
expanded as a Laurent series including double poles (~ 1/e 2 ), single poles (~ 1/e), and finite terms 
(~ e°). Further terms of 0(e) are unimportant here and will be left out. 

All results for the final-final and for all other dipole types can be found in [T]. 
2.5.2 Full result 

Having at hand the integrals for each dipole function, all individual dipoles present in a specific 
process can be collected to yield the overall infrared divergence of the subtraction term. Then, the 
starting point for the calculation of jet cross sections in the dipole subtraction formalism reads 



NLO 



E 

{m+l} 



m+l 



& a {m+l} | e=0 



da 



{m+l} |e=0 



+ 



/ 

J m 



{m} 



{m+l} 



T {m+l} 



,(47) 



e=0 



where ^{m+i} denotes the sum over all parton-level processes. However, the important point here is 
to exactly cancel the poles of the corresponding individual one-loop parton-level processes, which is 
done exclusively for each momentum and flavour constellation. Therefore, for each specific m-parton 
process at NLO only a selection of dipole functions related to (m + l)-parton processes contributes 
to the cancellation of the virtual divergences. In jl] it has been shown that this amounts to an 
effective reordering of phase space integrals and sums over parton configurations, such that 



a 



NLO 



E 

{m+l} 



m+l 



aa {m+l} |e=0 



{m} 



da 



{m} 



da 



{m} 



e=0 



(48) 



where da^ m ^ is the integrated dipole term that collects the integrals of all dipole functions and thus 
cancels the singularities of daY x . It is explicitly given by 

1™; 



da 



{m} 



da 



B 

{m} 



(49) 



where da^ m ^ x 1(e) is a shorthand for the following procedure: Write down the expression for da 
and replace the corresponding squared Born-level matrix element 



B 



\M {m} \ 



a (l, . . . ,m|l, 



,m) 



(50) 



with 



m (l, . . . ,m|I(e)|l, . . . ,m) m , 
using the insertion operator 1(e) as defined below. 



(51) 
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Finally, the full result for the integrated dipole term and the collinear counterterm as defined in Eq. 
(|14p for the most general case with hadronic initial states reads 

d<Tab(Pa,Pb) + da^,(p a ,p b ,n 2 F ) = [da^ b (p a ,p b ) x 1(e)] 



+E 



dx 



K a ' a ' (x) + P a ' a ' (xp a , x- fi F ) ) x da$ b (xp a ,p b ) 



+ J2 dx (K b > b \x)+P b > b '(xp h ,x^ F ))xdaZ,(p a ,xp b ) 



(52) 

where a and b again specify the initial state partons. The summation over a' and b' runs over all 
parton flavours, i.e. it includes gluons, quarks and anti-quarks occurring in the PDF. 
The insertion operator I reads 



(53) 



where the indices I and J run over initial and final state partons. The universal singular functions 
V/(e) depend merely on the flavour of I and are given by 



V,(e) 



V fl (e) 



C F 
C A 



2 ! ( ^CU - ^KiV; ) - + C A 



/50 _ 7f2 

I 9 2 



\ -T R N f ^- + 0(e) 



(54) 



6 3 '/ e 

with Aj being the number of contributing quark flavours. 

The complete singular structure in Eq. (i52j) is contained in [da^ b (p a ,p b )xl(e)} and the sum [da^ b (p a ,p b )x 
1(e)] + da^ b (p a ,p b ) must be finite for e — ► 0. 

The finite insertion operators K and P are given by 



K a ' a (x) 



as 
2vr 



i? aa (x) -ifp a s .(x) 



+d° a X T *- T «§ 



1 — x 



5(1 - x) 



T/, • T„/ 



T 2 

a 



-if a ' a (x) 



(55) 



and 



-'(w w ^^(4et,.t,,4 



(56) 



Note that here the index i runs over final state partons only. The flavour-dependent functions 

R aa R a,a ^ &nd pa,a ^ 

are defined in Appendix El As already mentioned, the factorisation- 
scheme dependent function Kp a s (x) vanishes in the commonly used MS-scheme. 

To obtain the final result for processes with no initial state partons only the I-term needs to be 
considered in Eq. (|52p . For processes with one initial state parton only, the result is obtained by 
using the I-term and one of the two integrals over K and P only, while omitting the contribution of 

K a ' a '(x). 
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2.6 Freedom in the definition of dipole terms 



As stressed before, the singular limits of the dipole functions are fixed by the requirement to cancel 
the singularities of the real correction matrix element. However, away from this limit there is some 
freedom for modifications. 

One possible modification has been presented in [21] . where a parameter a has been introduced 
which cuts off a dipole function for phase space regions far enough away from the corresponding 
singularity. Its main advantage lies in a significant reduction of the average number of dipoles terms 
to be calculated for each phase space point of the (m + l)-parton phase space of the real correction 
term. This constitutes an important alleviation of the calculational burden, since the total number 
of dipole terms grows approximately as m 3 . The a-modified subtraction terms also allow nontrivial 
checks of the implementation, since the total result must be independent of a. 

The a-modified dipole functions have been defined as follows: 

T> ij,k 0(a ~ Vij,k) , 
Vfj 0(a - 1 + x ij>a ) , 

Vf 6(a - Ui ) , 

V ai ' h 6{a - Vi) . (57) 

They will be employed later, in the implementation presented in this paper. Of course, such a 
redefinition of the splitting kernels also requires a recalculation of their integrals. The new a- 
dependent insertion operators I and K have been presented in |21j . 

Another simple modification is the addition of finite terms to the splitting functions, such as 

Vij,k + Vij,k * C , 
V§+(1- X ij:a ) * C , 

V k ai + Ui *C, 

V ai ' b + Vi*C . (58) 

The constant C directly ends up as a finite term in the integral of the splitting function and thus 
it can be easily included in the insertion operators of I and K, too. This again allows checks of the 
implementation, but it can also be employed to improve the numerical behaviour of the phase space 
integrals and to reduce the number of negative events. 



TV 
V' a 

v ,ai 



V' a 



v 



laiJ 



3 Implementation in AMEGIC++ 

The Catani-Seymour dipole subtraction terms have been implemented in full generality into the 
automatic matrix element generator AMEGIC++, based on it's version 2.0 |52pl . In particular this 
translates into AMEGIC++ being able to automatically generate all relevant parts of the NLO matrix 
element within the subtraction method except for the virtual matrix element. It can be applied to 
any process with massless partons for which the real correction ME can be generated, an extension 
to allow also for massive particles is foreseen. This includes standard model processes as well as 
implemented extensions, as long as there are no new strongly interacting particles involved. For 

3 A brief description of AMEGIC++ within the SHERPA framework can be found in |53| . whereas a full documentation 
of the (partly obsolete) version 1.0 is given in [2] with some extensions and results discussed in [51]- [56| . An update 
on the helicity formalism as it is used in the current version is documented in [54] . 
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standard model processes the boundary is currently at about six-eight partons (initial and final 
state). 

The new implementation aimed at a maximal reuse of already developed automated methods of 
amplitude generation and process management. Therefore, first a brief overview over the relevant 
parts of the code are given before the new implementation is described in some detail. 

3.1 ME generation 
3.1.1 Amplitudes 

The matrix element evaluation in the C + +- code AMEGIC++ is based on the evaluation of Feynman 
amplitudes using a helicity method based on the developments in [57]- [59]. The fundamental 
idea of this method is to introduce a helicity basis, in terms of which all physical spinors can 
be expressed. This allows to compute each amplitude directly as a complex function of physical 
momenta and helicity /spin states instead of computing traces of spinor products and 7-matrices for 
squared amplitudes. The colour within any amplitude is treated separately, i.e. in the first step all 
colour factors (the SU(3) structure constants f abc and i?-) corresponding to the fc-th amplitude 
are collected in an array 

Squared matrix elements can thus be written as 



is generated once for each process. In an initialisation run, AMEGIC++ generates all formulae for 
the amplitude calculation of the user-specified parton level processes including the colour matrix c%j 
(using a set of replacement rules for the colour algebra), it simplifies the expressions for the helicity 
amplitudes by identifying and factoring out common pieces and finally stores everything in C++ 
libraries. 

3.1.2 Process structure and organisation 

Since typically many parton level processes contribute to jet cross section calculations, usually a 
long list of processes needs to be computed. The corresponding structure in AMEGIC++ is as follows: 

• Any parton level process is represented by the class Single_Process, 

• Process_Group contains a (possibly recursive) list of such processes or groups of processes. 

All parton level processes sharing a specific common set of properties are grouped together in two or 
three levels of groups. In many cases there are subprocesses contributing to the same jet cross section 
which are very similar. Therefore AMEGIC++ applies a procedure to identify such processes in order 
to save computer resources and accelerate the calculation. The following checks are performed: 

• Direct comparison of amplitudes: check for processes that have identical graphs, where all 
involved particles have the same masses, widths and underly the same interactions (with cou- 
pling constants that differ at most in a constant factor). 

Example: QCD processes that differ in quark flavours only. 




(59) 



and hence a colour matrix of complex numbers 




(60) 
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• Numerical comparisons: check if the numerical result for a squared matrix element at a given 
phase space point is the same. 

Example: a quark is replaced by an anti-quark w.r.t. to the other process. 

For a set of processes that can be identified by this it is enough to compute one to know them all. 
In such a case, the corresponding matrix element squared is calculated only once and then recycled 
by the other processes. 



3.2 Generation of CS dipole terms 
3.2.1 Colour and spin correlations 

The starting point of the Catani-Seymour algorithm is detailed in Eqs. (148)) and (|49l) . supplemented 
with expressions like the one in Eq. (|22p for the individual dipole subtraction terms. The latter 
states that for any given process the Catani-Seymour dipole subtraction term for the real {m + 1)- 
parton correction term consists of the corresponding m-parton matrix element at Born level plus an 
additional operator that acts on colour and spin space. For the latter, only the limit e — ► needs to 
be considered. 



Colour operator: 

In all four dipoles, Eq. ([25]) . (f3Tj) . ([35]) . and (|39|) colour-correlated tree-amplitudes of the form 
\M^\ 2 = ^l,...,m|TVT fc |l,...,m) m (61) 



occur, where i labels the emitter and k the spectator. Denoting the colour indices of the 
external legs of the tree process explicitly by a,j and bj, this can be cast into 

= ^ ...i a > ...k a " ...m a ™\5 aibl ...T c aihi 

where T£ b = if acb , if the associated particle is a gluon, and T?- = t?-, if the associated particle 
is a quark. In other words, the colour structure for dipole terms can be generated by adding a 
gluon connecting the emitter with the spectator as illustrated in Fig. [2J The colour matrix for a 
dipole term is recomputed after this insertion using the available evaluation tool in AMEGIC++. 

• Spin space: 

For a quark splitting all spin-matrices are just proportional to 5 SS /, translating the quark spin 
to be exactly the same as for the Born- level m-parton matrix element. 

For the case of a gluon splitting, however, there are non-trivial correlation matrices. All of 
them can be cast into the generic form 

V^ = (fi\V\u) OC -gK» + ?¥L (63 ) 

where B and p are functions of the kinematic variables and momenta of the corresponding 
splitting. Their values are listed in Table HJ 



...T: kbk ...5 amb Jl bl ---i b '---k bk ...m b ™) m , 

(62) 
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1 




m+1 



Figure 2: Construction of the colour matrix for dipole terms: a gluon connects emitter and the 
spectator. 



dipole type splitting: B 



5 — 99 ~v _ ~* 1/(4^-) 

Z IPI z jpj , 

9^99 ^- i-jjd-^o - i-^d 1 -^) ,) /(2^\ 



FF 



FI 



5 — 99 ~u _ 9 * V(4£*j) 

gg ( 2 ~ i-jj+(i-gij,a) ~ l-gj+ci-aij,,)) /( 2 ^j) 



IF 



TT 5^99 vf-v-lf —Ab/( l - X i,ab) 



g — gg 



Table 1: Values for the functions defined in Eq. (j63f) . The variables are defined in the corresponding 
sections I2.4.1||2~T^1 The dipole type FF refers to the case where emitter and spectator are final 
state partons, IF refers to the case where the emitter is an initial state parton and the spectator a 
final state parton, etc.. 
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The structure of the splitting tensor as given in Eq. (j63[) is very similar to the polarisation sum for 
massive vector bosons in unitary gauge, except for the factor B and the fact that p can be timelike 
or spacelike. This analogy can be used to replace the tensor by a polarisation sum, i.e. 



B p A 



(64) 



Here the summation index A runs over four values, +, — , I and s. £ A is a sign that cannot be 
absorbed into the polarisation vectors e^. For a gauge boson with momentum 

p 1 = (j)o, \p\ sin#cos</>, \p\ sin sin 0, \p\ cosO^j , (65) 

the polarisation vectors are defined as 



1 



V2 



p2 



(0, cos 9 cos (f>Ti sin cf), cos 6 sin <fi±i cos <p, — sin I 



\f>\iPo-^: 



1 — B 



and the sign factors are given by 



i, e 



-1 if p 2 < 
+1 if p 2 > 



'+1 if p 2 < 
-1 if p 2 > ; B > 1 
+1 if p 2 > ; B < 1 



(66) 



(67) 



In order to calculate the dipole matrix element, the polarisation vectors of the splitting gluon are 
then replaced by the ones defined above. 



3.2.2 Organisation and process management 

To construct all dipole functions necessary to cancel the infrared divergencies of a given parton level 
real-correction process firstly all pairs of partons have to be determined that might emerge from the 
splitting of an emitter parton (initial state partons are charge conjugated for this procedure). This 
might be any quark (or anti-quark) and a gluon, two gluons or a quark and an anti-quark of the 
same flavour. Secondly, each of those pairs is combined with any possible third parton (acting as 
spectator) to define all possible dipole functions. 

Any individual dipole function is thus specified by: 

1. type (the specific combination of initial and final state for emitter and spectator), 

2. the specific flavours involved in the splitting, and 

3. the corresponding m-parton matrix element and its emitter and spectator particles. 
In order to construct the individual dipole functions, given by 

V = AiCljA^Fi...) . (68) 
the following ingredients are necessary: 
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1. A rule to map the (m + l)-parton phase space onto an m-parton phase space. 

2. The corresponding splitting function for the dipole. This consists of two parts, a scalar function 
F{...) of the kinematic variables of the splitting and a spin correlation matrix. As discussed 
above, for quark splittings the matrix is simply <5 SS /, for gluon splitting the matrix is represented 
by an outer product of pseudo-polarisation vectors, which are also functions of the kinematic 
variables of the splitting. 

3. The colour matrix C^-, respecting the extra colour correlation. 

4. Amplitudes A, of the corresponding m-parton matrix elements. For gluon splitting cases these 
amplitudes have to be calculated replacing polarisation vectors of the splitting gluon by the 
pseudo-polarisation vectors introduced above. 

The calculation of any dipole function is organized in the class Single_DipoleTerm, each instance of 
this class representing one dipole. This class controls the ingredients for the calculation: Firstly there 
is a Born-level m-parton matrix element of the original AMEGIC++ implementation, just extended 
such that it includes the additional colour correlation. Secondly there is a class Dipole_Splitting_Base 
that completely organizes the splitting function itself. Specified by the type of the dipole (initial and 
final states for emitter and spectator) and the type of the splitting (determined by the contributing 
flavours) it takes care of the mapping between the m + 1-parton and the m-parton phase spaces 
and of the calculation of the splitting function (including the polarisation vectors to encode the spin 
correlation) . 

Above that the class Single_Real_Correction handles all contributions to an infrared regular- 
ized parton level process. This consists firstly of an (m + l)-parton tree level matrix element in 
the original AMEGIC++ implementation. Secondly it contains a list of single dipole functions, sim- 
ply determined by looping over all partons and selecting valid dipole configurations. The classes 
Single_Real_Correction and Single_Process are derived from a common base class in a way such 
that the class Process_Group can be reused to also organize the infrared regularized parton level 
process in groups of common features up to all subprocesses contributing to a jet cross section. 

Similarly to the case of tree level processes in AMEGIC++, also here a mapping of parton level 
processes that lead to identical or proportional results can be used to speed up the calculation 
and save computer resources. To this end, the following automatic identification strategies are 
implemented: 

• If two real correction processes can be mapped (using strategies described in section I3.1.2j) 
then also the whole Single_Real_Correction is mapped. 

• For single dipole terms a unique identification algorithm proceeds as follows: Two terms can 
be mapped if the included m-parton process can be mapped and if the three particle labels 
(numbering the the external particles of the real correction process) to identify a dipole are 
identical. 

• Many of the born matrix elements within the dipole terms will be identical. However, since 
different dipoles require different momentum mappings they have to be recalculated. Only the 
calculation routine can be shared. 
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3.3 Generation of the finite part of integrated dipole terms 



3.3.1 Analytical structure of the full result 



The starting point of the discussion of the finite pieces of the integrated dipole terms is Eq. (I52|) . 
where now the phase space integration as well as the summation and integration over the incoming 
parton flavours and momenta is made explicit. Then, terms inside the m-parton integral come 
from subtraction terms integrated over the phase space of the extra parton emission and from the 
collinear counterterm for the general case of a NLO cross section with initial state partons. The terms 
inside the (m + l)-parton phase space integral in contrast corresponds to the dipole subtraction bit. 
Altogether, and including the convolution with parton distribution, the relevant term to be evaluated 
can thus be cast into 

^2 / d^d^ f aim ,^f) Mm ,^ 2 f)\ / &°ab{mp,mP) + / ^abimp^mP^V, 

a ^ J l J m+l J m 



a.b 



/ dmdmfa(vi,VF)fb(v2,Li F ) \ / [dof fe (f?ip, mv) x 1(e)] 



dx 



K a ' a ' (x) + P a ' a ' (xmP, x; /4) ) x da$ b (xrnp, mP) 



+ J2 dx \(K b >"(x)+P b '"(xT}2p,xitJ? F j) xdag,(mP,xmP) 



(69) 

The only correlation of the insertion operators I, P, and K with the Born level matrix element is 
within color space. To be more specific, this implies that only the following structures emerge 



da 



d(T ab (Pa,Pb) 
(Pa,Pb) 



B(i,j) 
ab 



t (l, . . . ,m;a,b\\l, . . . ,m;a,b) m and 
t (l, . . . ,m;a,b\Ti ■ Tj\l, . . . ,m;a,b), 



(70) 



for all i j, where i and j may label both final and initial state partons. Since any of the appearing 
matrix elements with insertion operators can be written as a sum of such structures, the color factors 
will be skipped in the following and the operators will be treated simply as scalar functions. 

The terms P and K induce dependences on x, which combined yield result in the structure 



(g(x))+ + 5(1 - x)h(x) + k(x) . 



(71) 



Here, h(x) and k(x) are regular functions in x and the '+ '-distribution is defined by its action on a 
generic test function a(x) 



I dx a(x) (g(x)) + = / dx [a(x) — a(l)] g(x) 
Jo Jo 



(72) 
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Then the r.h.s. of Eq. (J69J) can be cast into the form 
/ ^m^fafo, fJ? F )fbfo, Vf) 



+ V f dx (g a > a '{x) [da^xrup, m p) - da$ b fo p , m p)] + k a ' a ' (x)da^ b {xrnp, m p) 



+h a > a '(l)d<jZ b ( m p, m p) 
i 



+ V / dx (g b ' b ' (x) [da^y fop, xrj 2 p) - dcrf b , fop, mp)] + k b ' b ' (x)da B b , fop, xrj 2 p) 



+h b ' b (l)da^fop, m p) 



(73) 



The functions g a ' a (x), k a ' a (x), and /i a - a (1) can be read off the corresponding functions in App. |A"1 

Computationally the most demanding part is the actual Born- level cross section der^, due to its 
potentially expensive multi-particle matrix element, which typically suffers from factorial growth 
with the number of external particles. Thus, the calculation can be significantly accelerated if the 
expression is rearranged such that da^ b has to be computed only once for a single configuration at 
a given phase space point. This can be achieved by changing the integration variables rj to rf = xrj. 
After renaming 77' back to rj and reordering the summation over a and a' (b and b') the expression 
above reads 



/ driidri 2 f a fo, (j, 2 F )f b fo, (j? F ) / da^fop, rj 2 p) x < 1(e) 

a,b J Jm { 

+ £ /' dx [M2^M (/, >} + ,.> (x) ) _ M2i41/.« (I) 

LitMH,/*?) V J fafo,V F ) 

~ fafo,H F ) V 7 

+ ^/>[^f him, A) 



(74) 



where the G a,b (i]) = dxg a,b (x) are analytically computed. 

The insertion operator 1(e), Eq. (|53p is given as a Laurent series in e. For the implementation the 
interesting part is oc e°, since the poles must have been analytically extracted befor^l- 



4 For testing purposes, however, it is trivial to also determine the coefficients of the e 2 - and e 1 -poles and to 
compare with known results of virtual correction terms. 
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3.3.2 Implementation and Organisation 

The numerical calculation of the finite contributions from integrated counterterms is organized as 
Eq. (174)) suggests, i.e. the basic unit (class Single_Virtual_Correction) covers everything that is 
associated with a specific m-parton cross section. 

For the actual calculation, basically all colour correlated matrix elements in Eq. (|7Up are necessary. 
The contributing amplitudes are, of course, the same for all of them, only the colour matrix is 
different. Therefore, a generalized version of Single_Process is employed that is able to deal with 
a multitude of colour matrices to calculate all required matrix elements at once. Anything else 
needed for the calculation of the finite contribution is a long list of rather simple scalar functions 
and constants. The integration over x is done numerically, i.e. for each set of external momenta x is 
diced within the corresponding interval. 



3.4 Phase space integration 

Together with the automatic generation of matrix elements AMEGIC++ also generates specific, 
process-dependent phase-space mappings for efficient integration. There, some a priori knowledge 
about the integrand is used [60]- [62] together with self-adaptive Monte Carlo integration meth- 
ods |16|63] - |65j . Here, the general method will briefly be summarized for the case of tree-level 
processes and its application to the integrals coming with the subtraction method will be discussed. 



3.4.1 Importance sampling and Multi-channel integration 

The general idea behind importance sampling is to improve the numerical behaviour of an integrand 
by a change of integration variables, 

/ f{x)dx = / dy , where - = — — . (75) 

J J g{x{y)) 9 dy 



The new variable y is chosen in a way such that £ is a sufficiently smooth function, leading to 



I 

9 

a reduced error estimate of the integration. Typically, the weight g is chosen as a simplifica- 
tion/approximation of /, such that the integral y = J gdx can be analytically solved. 

For phase space integrals maps relating vectors of uniformly distributed random numbers {aj} inside 
the interval [0, 1] to the four-momenta of the external particles of a physical process {pj} 

{ Pj } = X({ ai }) (76) 

are in the center of the sampling process. The weight function in such a case g is then determined 
by 

1 d* n (X ({<*})) 

9 d{a t } ' 1 > 

For a single squared amplitude it is easy to determine suitable momentum mappings and weights: 
invariant masses of the propagators are determined according to the propagator and angles for 
particle splitting are chosen isotropically; the finial state momenta are then determined out of those 
variables. For instance, for a massless propagator the invariant s would be generated by 

s = K^ + Cl-a)^] 1 ^ , (78) 
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with the corresponding weight 



9 = (79) 

The constants s max and s m i n are upper and lower boundaries of the invariant mass, which depend 
on the overall topology of the phase space point and potential cuts; v in contrast is an effective 
exponent for the propagator, subject to choice. 

Weight distributions for contributions from several amplitudes can be then combined using the 
multi-channel method. A total weight function G is defined through 

G = Y, a ^9k, (80) 

k 

where the g k are the weight functions for the single contributions (channels) and the a k are arbitrary 
coefficients with a k > and Y^k a k = ^he corresponding momentum mapping is then given by 

fc-1 k 

X({a l },o;) = X k ({a,i}) , for ^ a/ < a < ^ a/ . (81) 

1=1 l=i 

The multi-channel method relies on automatically adapting the coefficients a k such that the variance 
of the phase space integral is minimized. 

Further refinement 

The efficiency of the integrator is improved if additionally the self-adaptive VEGAS algorithm [61] is 
applied on the channels. VEGAS is very efficient in the numerical adaptation to functions, where the 
peaking behaviour is not too extreme and which are factorisable to a product of one-dimensional 
functions. This is clearly not given for full matrix elements. However, the structure represented by 
a single channel fulfills this condition. Thus, in AMEGIC++ VEGAS is used to adapt selected channels 
to structures that go beyond their rough approximations and which are typically hard to specify 
analytically or which are a priori unknown. 

For each channel VEGAS is used to generate a mapping £ from uniformly distributed random numbers 
to a non-uniform distribution, still inside the interval [0,1], and a corresponding weight v^. To 
combine this with the multi-channel method the mapping X({ai}) for single channels must meet 
the requirement to be invertible. The full map reads 

fc-i k 

X({fli},a) = X k (£ k ({ai})) , for < a <^ai . (82) 

1=1 i=i 

For a momentum configuration {pj} the weight is therefore given by 

G({Pj}) = Y, a ^k{{Pj})v k {X-\{ Pj })). (83) 

k 

Subtracted processes 

The subtraction method necessitates the evaluation of two independent integrals, namely integrals 
over the m-parton and the (m + l)-parton phase space. In both cases mappings generated for the 
tree level process of the same dimensionality are used. 

For the integration of the (m+ l)-parton phase space soft and collinear regions must be included. In 
this case the lower limit for the invariant masses of many propagators (e.g. Eq. (|78|) ) must be zero. 
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To keep the integral over the weight finite the exponent v must be set to a number smaller than 1. 
The actual shape of those propagators is hard to specify a priori. It depends on the jet definition 
and on the balance between the real correction process and the subtraction term (the integrand can 
be positive or negative). Taken together, however, it seems not unreasonable to assume a small 
exponent. The VEGAS refinement adapts very good to the actual shape and the final integration 
efficiency after optimisation has only a weak dependence on the initial value of v. Since the VEGAS 
algorithm optimizes on the variance of the integrand it can, to some extend, also deal with the 
numerical problems related to "missed binning" , which will be discussed in the following section. 

The m-parton phase space is much simpler. Since most parts of the integrand are proportional to 
the born matrix element it tends to work very well with this phase space setup. 

3.5 Cuts and analysis framework for NLO calculations 

Triggers and observables for NLO calculations have to be chosen with care. The general strict 
requirement not to spoil the cancelation of infrared divergencies has already been discussed in section 

E3 

Before going into any details concerning cuts, it is important to notice that a rule is mandatory 
of how cuts act on the different contributions to the NLO cross section. This rule must exist in a 
m-parton and a (m + l)-parton version, where the latter needs to satisfy the conditions of infrared 
safety in degenerate phase space regions. In practical terms, this implies that the (m + l)-parton 
version of the cut must allow for exactly one parton to become soft or collinear, while the m-parton 
version has to omit all singular regions. 

Second, Eq. ([TH]) requires for the cut of the m+1 phase space integral to be applied separately to the 
real correction process (using the m + 1-parton version) and to each dipole term (using the m-parton 
version, applied on the momenta of the mapped m-parton configuration). In general there might 
be kinematic configurations, where the real correction process ends up outside the accepted phase 
space region but some dipole terms do not and vice versa. This leads to the problem of "missed 
binning": if such a configuration occurs close to a singular region, large contributions result, which 
do not cancel completely. Ultimately, this leads to large numerical fluctuations, which need to be 
addressed. This is a common issue for all subtraction methods. 

So far, the following cuts have been made available in AMEGIC++: 

• A simple cut for jets is implemented as follows: a suitable jet algorithm (e.g. kx) [HS]- [69J is 
used to construct jets from the final state partons and their momenta. Then the number of 
jets above a given py-cut is counted. A phase space point is valid if this number is greater or 
equal m. 

• Of course also cuts that only act on particles not taking part in strong interactions can be 
applied. If initial-initial dipoles are present this also has to be done separately for the real 
correction and for the dipole terms, since the momentum mapping in this case modifies all 
final state particle momenta. Implemented are cuts on invariant masses, on total or transverse 
energies, on rapidities or on particle angles w.r.t. the beam. 

Sherpa's ANALYSIS-package has been extended to be able to deal with weighted events from the 
NLO subtraction procedure. For example, and to be more specific, consider the case of a cross 
section which is differential to some infrared safe quantity F, i.e. a distribution to be binned in a 
histogram d-F. For the m-parton integral no special treatment is mandatory: for a given momentum 
configuration, dF can directly be evaluated and filled into the corresponding bin. For the real 
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correction and the dipole subtraction functions in the (m + l)-parton integral, F has to be evaluated 
for each contribution separately, similar to the phase-space cut. Again, the problem of "missed 
binnings" appears, if contributions to a single event end up in more than one bin. 

4 Checks of the implementation 

In this section a number of tests of the correct implementation of the subtraction algorithm and of 
the integration routines are described. These tests are mainly technical in nature, results relating 
to truly physical observables are discussed in the next section, Sec. EJ 

4.1 Explicit comparisons 

Before moving on to technical checks, it is worth stating that a number of direct comparisons of 
individual terms from the program presented here with those obtained from M. Seymour's Fortran 
code DlSENT have been performed. The latter is a dedicated program to compute NLO cross sections 
for the deep inelastic scattering processes e~p — ► e~ +jet, e~p — ► e~ + 2 jets and for electron-positron 
annihilation to two and three jets. This direct comparison is possible, since DlSENT uses exactly the 
same subtraction formalism, allowing to compare individual terms at given phase space points. All 
terms listed in the following showed full agreement of the two codes, up to the numerical precision. 

The comparison included: 

• Dipole subtraction terms for the real correction: 

all flavour configurations for dipoles with final state emitters and spectators as well as for 
dipoles with initial state emitters / final state spectators and final state emitters / initial state 
spectators have been checked. 

• Terms from the finite part of the insertion operator I, 
cf. Eqs. (gS]) and (|53|). 

• Terms from the insertion operators K and P for the case of one initial state parton, 
cf. Eq. (|52|) and the implemented version Eq. (|74|) . 

Furthermore integrated results of the virtual and real parts of the NLO corrections in this subtraction 
scheme where compared and agreed within statistical errors for all accessible processes. 

4.2 Test of convergence for the real ME 

An obvious first technical check of the overall package consists of testing the convergence behaviour 
of the dipole subtraction terms close to the singular region. To this end, the m + 1-parton phase 
space of the regularized real correction part is numerically integrated over. The crucial issue is to 
ensure that the integrand remains finite over the full phase space, in addition the performance of 
the integration algorithms deserve some consideration. 

Clearly, for the numerical calculation a small phase space region around each singular configuration 
has to be cut out. Although the dipole terms are expected to become equal to the matrix elements 
there, technically speaking infinite or very large numbers must be subtracted in this region, leading to 
large fluctuations and hence to errors due to the limited numerical precision at which the calculation 
is performed. Therefore a variable a m i n is introduced, which on the basis of kinematic variables of 



28 




-1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 -1 -2 -3 -4 -5 -6 -7 -8 -9 -10 -11 



log(a ut ) log(a ut ) 

Figure 3: Dependence of the subtracted real emission cross section on a cut for (a): e~e + — > 2 jets; 
(b): e~e + — > 3jets, both at a CM energy of 100 GeV; (c): e~p — > e~ + jet with a 50 GeV electron 
beam and protons at 500 GeV; (d): pp — ► 2jets at a CM energy of 14 TeV. To obtain a well-defined 
LO cross section for (b) at least two jets with a fc^ ur ' > 10 GeV, for (c) a transverse energy of the 
scattered e~ > 10 GeV and for (d) at least two jets with p± > 40 GeV are required. 
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Figure 4: Normalized absolute values of cross sections in bins of a m \ n . Setups and phase space cuts 
are the same as in Fig. [3l 
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corresponding dipole functions, reads as follows: 



a. 



min (adipole) 

dipoles 



(84) 



where 



Q-dipolc 



Uij,k 



Ui 
Vi 



for T>i jik 

for P« - 

for Vf - 
f or _ 



—dipoles (final state emitter, final state spectator) 
-dipoles (final state emitter, initial state spectator) 
-dipoles (initial state emitter, final state spectator) 
-dipoles (initial state emitter, initial state spectator) 



.(85) 



This parameter a serves as a cut-off in such a way that for an externally given parameter a cut 
kinematic configurations with a m i n < a cu t are omitted. 

In Fig. [3] the dependence of the subtracted cross section on a cut for four sets of real correction 
processes, namely e~e + — > 3jets, e~e + — ► 4jeis, e~p — ► e~ + 2jefs and — > Sjets. All types of 
dipoles and splitting functions contribute to the dipole terms which are necessary to regularize those 
processes. It is apparent that for a cut ~ 10 -5 the cross section stabilizes close to its final value. 

To study the numerical behaviour near the singularity in more detail, in Fig. 0] the absolute value 
of the subtracted cross section, binned in intervals of ctmin is depicted. For all studied processes the 
contribution to the cross section drops down by at least four orders of magnitude with decreasing 
a m i n and confirms the observations for full subtracted cross sections made before. 

The strong increase accompanied with statistical errors of 100% or larger for a m i n values below 
10~ 9 — 10 -11 signals defects due to the limited numerical precision (double precision ~ 10~ 12 ). One 
reason is the already mentioned numerical problem when subtracting extreme large and almost equal 
numbers. Another reason is the precision of the momentum four-vectors itself, because the precision 
of the external particles residing on their m = mass shell is also limited by the numerical precision. 
This of course may consequently lead to errors of that order in the matrix element calculation. Thus, 
Fig. 2] allows to determine best choices for a cu t, somewhere between 10 -9 and 10~ . 



4.3 Consistency checks with free parameters 

In section [2~B1 ways of modifying the subtraction terms without changing the singular behaviour have 
been discussed. Such modifications can be employed for non-trivial tests of the implementation, since 
the modifications will affect both, the real part and the virtual part of the NLO cross sections, with 
their sum remaining constant. 

In Fig. [5] the total NLO correction for the cross sections of e~e + — > 2jets, e~e + — > 3jets, e~p — > 
e~ + jet and pp — > W~ — > e~v e and their real and virtual contributions are displayed as functions 
of the parameter a, as introduced in section 12.61 The fact that the sum remains constant within 
statistical errors provides a non-trivial confirmation of the correct implementation of the algorithm. 
It should be noted here that the calculation of the cross section of the processes under consideration 
invokes all types of dipole functions as well as the most general case of the insertion operators from 
the integrated dipole terms. 

By using the same number of phase space points for each integral and comparing statistical errors, it 
can be seen that this parameter can also be used to optimize the numerical behaviour. Clearly, best 
results are obtained if the values of the virtual and real contributions are both as small as possible, 
thus reducing the size of the fluctuations. It should be noted that the error bars in Fig. [5] are given 
not including the leading order part of the cross section. Since relative errors for the latter can be 
expected to be much smaller if evaluated for the same number of phase space points, the (relative) 
statistical error for the full NLO cross section will be significantly reduced. 
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Figure 5: NLO corrections as a function of the parameter a in the definition of the subtraction terms 
(see Eq. ([FT]) ) for the total cross sections of (a) e~e + — > 2jets, (b) e~e + — ► 3jets, (c) e~p — ► e~ + j'ei 
and (d) — > W~ — ► e^P e . The results and error bars in the difference plots are determined after 
calculating 500000 phase space points of the real contribution, which typically dominates the total 
statistical error. 
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5 First physical applications 



In this section, some simple applications will demonstrate the performance of the dipole subtraction 
procedure, as implemented, for the calculation of physically relevant observables. The born matrix 
elements, dipole subtraction terms to regularize the real correction and corresponding finite terms 
to be added to the virtual correction were generated automatically by AMEGIC++. The one- loop 
amplitudes have been explicitly implemented for the considered processes. 

5.1 Three-jet observables at LEP 

To compute three jet cross section at next-to- leading order the one-loop matrix element given in 
|36j has been implemented. The expression given there is averaged over the direction of incoming 
momenta, which is sufficient for observables that are not correlated to the beam direction. 

In Fig. [6]LO and NLO predictions are displayed for observables sensitive to 0{as)- In particular, 
the event shape observables 1-Thrust, Major, C-parameter and the Durham 3^2 jet rate are 
compared with measurements performed at LEP on the Z°-peak by DELPHI [70]. All data are 
normalized to unity. The normalisation for the calculated cross section, however, is somewhat 
complicated. This is because in the calculation three-jet events are required in each case, translating 
into the necessity to apply a phase space cut. On the other hand, the data are more inclusive and 
also include comparably soft regions, where fixed-order perturbation theory is known to fail and 
must be supported by resummation techniques. The normalisation for the calculations has thus be 
chosen such that it agrees with data in the "safe" regions. This exposes the differences between 
LO and NLO calculations best. As a consequence, the corresponding normalisation factor of both 
calculations is not identical. From the results of Fig. [6] it can be deduced that for all observables the 
range described sufficiently well by the calculation is extended for the NLO calculations. For both, 
the region described by soft physics (left side in all plots) as well as phase space regions populated by 
additional hard QCD radiation (right side in all event shape plots) the prediction has been improved. 

5.2 DIS: e-p -> e~ + jet 

The one loop matrix element for this process is given by the well known expression 

lMiU v) = IMI^^^ (^)*{-|-|-8 + 0(e)) , (86) 

where Q 2 = —q 2 > with q the momentum transfer between the electron and the proton. 

Fig. [7] shows differential cross sections w.r.t. the transverse momentum and rapidity of the scattered 
electron and the hardest jet at leading and at next-to-leading order. The CM-energy has been taken 
as V 10 5 GeV, corresponding to a 50 GeV electron beam and a 500 GeV proton beam. The parton 
distribution function CTEQ6M |71j has been employed, factorisation and renormalisation scales 
have been fixed to Q 2 . A phase space cut on the electron of px > WGeV has been imposed. The 
NLO correction for this setup is comparably small, for the total cross section it is of the order of 5% 
and negative. The ratios of NLO and LO calculation, however, are not constant for all observables. 
At NLO the cross section rises for increasing momentum transfer, up to a correction of 40% for 
transverse momenta of electron and jet of the order of 150 GeV. 
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Figure 6: The event shape observables 1-Thrust, Major, C-parameter and the Durham 3 — * 2 jet 
rate at LEP I compared to measurements by DELPHI [70]. The LO and the NLO predictions have 
been normalized to the data separately in a region, where agreement can be expected. The dashed, 
dashed-dotted and the dotted lines are the Born, real and virtual contribution to the NLO cross 
sections, respectively. 
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Figure 7: Distribution of transverse momentum (left plots) and rapidity, denned in the beam CM 
frame (right plots), of the scattered electron and of the hardest jet in deep inelastic scattering, 
calculated at leading order and next-to-leading order. The CM-energy has been taken as VW 5 GeV. 
A phase space cut on the electron (pr > WGeV) has been applied. For the rapidity distribution 
of the first jet a px > 15GeV has been required. Dashed and dotted lines denote the real and the 
virtual corrections to the Born cross section, respectively. The lower panels of each plot show the 
ratio between the leading order and the next-to-leading order results. 



35 




10 12 3 u -> u 1UU 

r((W") P x (e") [GeV] 



Figure 8: Rapidity distribution of the VF~-boson (left plot) and transverse momentum of the electron 
for the process pp — > W _ — ► e~v e at Tevatron Run II, calculated at leading order and next-to-leading 
order. Dashed and dotted lines denote the real and the virtual corrections to the Born cross section, 
respectively. The lower panels of each plot show the ratio between the leading order and the next- 
to-leading order results. 



5.3 W production at Tevatron 

The one-loop virtual contribution to this process can be obtain by crossing relations from Eq. ([86 
and is given by 



Ma-fad = \M\ {barn) 27r t{1 _ £) {-qt) \-^-- e - 8 + 7r +a{ ' 

where now Q 2 = s, the CM energy squared of the incoming partons. 

Fig. [8] shows cross sections for Tevatron Run II, differential in the rapidity of the VF~-boson and the 
transverse momentum of the electron, respectively. The parton distribution function CTEQ6M [7T] 
has been employed, factorisation and renormalisation scales have been fixed to rn^y. 

The total and differential cross sections are in full agreement with predictions obtained using the 
next-to- leading order parton level generator MCFM [20J. 



6 Conclusions and outlook 

In this publication a fully automated implementation of the Catani-Seymour dipole formalism in the 
framework of the matrix element generator AMEGIC++ has been presented. It allows to automatically 
generate the process-dependent real correction terms for given Born cross sections with massless 
external particles and the corresponding real subtraction terms. The integration of the subtracted 
real correction terms is performed automatically with a multi-channel method, giving rise to an 
appreciable convergence. The implementation has carefully been checked for correctness, invoking 
consistency checks with free finite terms which may be added to the subtraction terms. Through 
the explicit inclusion of virtual terms a parton-level calculator is so available. 
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In the future, the code will be further updated to include massive external particles and to provide 
a full parton- level generator at NLO. 
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A Insertion operators etc. 

In this appendix, the ingredients of the master equation Eq. (|52p . 
dVab(Pa,Pb) + da^ b (p a ,p b ,n 2 F ) = [da^ b (p a , p b ) x 1(e)] 

+ 12 [ dx [(K a ' a '{x) + P a ' a '{xp a ,x;LiF)) x da*, b (xp a ,p b ) 

a' J ° 

+ J2fdx [(K b ' b '(x) +P b > b '(xp b ,x;fi 2 F)) x da*,(p a ,xp b ) 



will be repeated, a and b specify initial state partons, and the sum runs over all accessible a' and b' 
occurring in the PDF. The insertion operator I is given by 



2 \ 6 

(89) 



cf. Eq. (|53p . and again the indices / and J run over all initial and final state partons, while the 
universal functions V/(e), encoding the singularity structure, merely depend on the flavour of / and 
read 

V/(e) = Tf (l-y)+7iQ + l)+tfj + 0(e), (90) 
cf. Eq. ([53]). The individual 77 and Ki will be listed in Eqs. ([Ml) and ([§5]) . 

The factorisation scale dependent terms are proportional to insertion operators P a ' a ' ({p}, xp a ,x; up), 
which read 

P a > a '({p},xp a ,x;vl) = ^p-'(x)-^ £ T/ .T a ,ln-^— . (91) 
The regularized Altarelli-Parisi kernels P ab (x) are listed in Eq. (|96p . 
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The factorisation-scheme dependent terms are proportional to the initial-state insertion operators 
K. For one initial-state hadron only, this operator reads 



K a > a '(x) = — { K aa \x) - K a F a ' s {x) + 5 aa ' V 



1 



1 — x 



+ 8(l-x) 



(92) 



with the functions ifp°g (x) and K aa> (x) given below, cf. Eqs. ()97[) and (|99p . and with the ji listed 
in Eq. (I94p . Note that the subscript "F.S." denotes the factorisation scheme. For two initial state 
partons, the initial-state insertion operator is given by Eq. (|55p . 



K a ' a (a) 



as 
2vr 



x)-K^{x) 



+d aa'^ T .. Ta ^ 



1 — X 



+ 6(l-x) 



T 2 

a 



(93) 



with the functions K aa ' (x) given in Eq. (|98|) . 

The 7/ and X/ occurring in Eqs. (|90p are related to integrals of the Altarelli-Parisi kernels listed 
below, Eq. (|96l) . and read 



3 n 11 r> 2 t at 

lq = lq = ifF , lg = ~^C A ~ ^T R N f 



and 



K q = K- q 



7 n u 



67 7T 2 



18 6 



10, 
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(94) 



(95) 



respectively. The Altarelli-Parisi kernels emerging in the factorisation-scale dependent terms of Eq. 
(j9H are 



P q9 {x) 
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(96) 
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The functions K ab (x) are explicitly given as 

K qq (x) = K qq (x) = 0, 
K q9 (x) = K q9 {x) 

K 9q (x) = K 9q (x) 
K qq {x) = K qq {x) 



P q9 (x)ln- — - + C F x , 

x 

ps9( x ) i n 1^1 + 2T R x{l - x) , 



C F 



x 

2 l-x 
In 



1 — x x 

2^ 



(1 + x)ln- + (l-x) 



X 



K"(x) = 



-<5(1 -x)(5-tt z )C f , 
2C A 



1 , 1 -x 
In 



1 — x 



, 1 — x (\ — X , ,_, . 
+ ln ( 1 + x(l - x) 



x 



-6(1 -x) 
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)C A -—T R N f 



(97) 



whereas the functions K ab (x) read 



K qq (x) = K qq (x) = , 
K q9 (x) = K q9 (x) = P q9 (x) ln(l - x) , 
K 9q (x) = K 9q (x) = P 9q (x) ln(l - x) , 

2 



l^(x) = K w (x) = C F 
K"(x) = C A 
+2 



1 — x 
2 



ln(l — x 
ln(l - x) 



7T 



<5(1 - x) - (1 + x) ln(l - x) 



7T 



l-x'" y ~ "'I , 3 
1 — X 



(5(1 -x) 



X 



— 1 + x(l — x)J ln(l — x) 



(98) 



Finally, the factorisation-scheme dependent terms are given through 



^DIS — "^DIS 



^DIS — ^DIS 

K 19 _ K qg 
-"■DIS ~~ ^DIS 

^DIS 

-"■DIS ~~ ^DIS 



= o, 

= C F 



l + x 2 

1 — X 



In 



1 — x 



9 + 5x 



(x 2 + (1 - x) 2 ) In i— ^ + 8x(l - x) - 1 



-"■dis ' 
. 



(99) 
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